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Abstract 

Classical wave equation is generalized for the case of viscoelastic materials obeying frac- 
tional Zener model instead of Hooke's law. Cauchy problem for such an equation is studied: 
existence and uniqueness of the fundamental solution is proven and solution is calculated. 
Ph ' 

ril ■ Mathematics Subject Classification (2000): Primary: 26A33; Secondary: 46F12, 74D05, 

; 74J05 

^ Keywords: fractional derivatives, Laplace and Fourier transforms, fundamental solution, 

■ fractional Zener model 

1 Introduction 

^ , The need for a use of derivatives of arbitrary real (or even complex) order, also called fractional 

' derivatives, is permanently underlined by various applied problems. Although, the idea about 

^>«^ , differentiating with respect to noninteger order originates from the work of Leibnitz and Euler at 

^SJ ' the turn of the 18th century, the last 30 years have brought a real expansion of fractional calculus. 

Certainly, it is due to great number of possibilities for its application in many diverse fields such as 
mechanics, viscoelasticity, automatic control, signal processing, stochastic and finance, biomedical 
engineering, etc. For a detailed exposition of the theory of fractional calculus we refer to [16] or 

m- 

In this paper we study existence and uniqueness of solutions for the system 



a{x,t)+ToDfa{x,t) = e{x,t) + oD^e{x,t), (1) 

, , d 
e[x,t) = —u{x,t), 

where a; G M, i > 0, cr, u and e are stress, displacement and strain, respectively, considered as func- 
tions of X and < t < 1 is a constant and qD", < q < 1, denotes the left Riemann-Liouville 
operator of fractional differentiation of order a. This system describes waves in viscoelastic media. 
System ([T]) is subjected to initial conditions 



d_ 

dt 

cr(x,0)=0, £(x,0)=0. 



u{x,0) = ua{x), —u{x,0)^vo{x), 
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as well as boundary conditions 



lim u{x,t)=0, lim a{x,t)—0. (3) 

a:— fitcjo X— ^±00 



As we shall show, system ([Ij can be reduced to the equation 



92 92 
— t) = Lit)—u{x, t), x e M, t > 0, (4) 



which will be called wave equation for fractional Zener type viscoelastic media or shortly FZWE 
(Fractional Zener Wave Equation) . Here L denotes a linear operator (of convolution type) acting 
on 5'(M), whose explicit form is given by 



L{t) = L 



-1 



1 + s" 
1 + rs" 



In fact, equation (|4]) will be the subject of our consideration. System ([T]), and in particular ([4]), 
generalize classical wave equation. 

Equation that describes waves occurring in the elastic media are obtained by using basic 
equations of elasticity (see [2])- Our interest is restricted to an infinite elastic rod (one dimensional 
body), positioned along x axis, that is not under influence of body forces. Equations of elasticity 
then read 

d 92 

— (7(a;,t) Pg^u{x,t) (5) 

a{x,t) = Ee{x,t) (6) 
d 

where a; € R, t > 0, u, u and e are stress, displacement and strain, as above, p = const, is the 
density of the media and E — const, is the Young modulus of elasticity. Equation ([5]) is the 
equilibrium equation and it is the consequence of the Second Newton Law. Equation ([6]) is the 
constitutive equation and in the case of elastic media it is known as the Hooke law. Equation Q 
is the strain measure for the local small deformations, i.e. connection between strain and displace- 
ment. Wave equation is now obtained by substituting equation ^ into ^ and subsequently (O 
into dS]). It reads 

Qf2^i^^t)-c^<x,t), c=y-. (8) 

Newly introduced parameter c can be physically interpreted as the wave speed. 

The first authors who introduced the fractional Zener model were M. Caputo and F. Mainardi 
[SI [?]• Such model was also considered in [T7] but without using fractional calculus. A similar 
problem of stress waves in a viscoelastic medium was investigated in [5] by using the inversion 
formula of the classical Laplace transform. Wave equation ([8]) can be generalized in the framework 
of fractional calculus in several ways. One of the possibilities is to substitute second order partial 
derivative with respect to time by fractional one, as it has been done in |13j for one dimension, 
and in [10] and [11] for d dimensions, € {1,2,3}. Time partial derivative can also be replaced 
with either two fractional derivatives of different order, as in [3], or by distributed order fractional 
derivative (integral over given range of order of fractional derivatives), as in [T^] and [13] for 
a e (0,1), and in (J and [5] for a G [0,2]. In each of these papers the intention was to stress 
similarity between wave and diffusion equation, since by allowing the order of fractional derivative 
to be a S [0, 2] both equations can be obtained as special cases. 

In this article wave equation ^ is generalized for the case of viscoelastic media described by 
fractional Zener model. Generalization is conducted so that constitutive equation ^ is changed, 
since Hooke's law is inappropriate to describe viscoelastic media. Other two equations, namely 
equation of motion of deformable body ([5]) and relation between strain and displacement ([7|), 
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are valid for any type of deformable body, with mentioned restrictions, therefore they remain 
unchanged. Similar generalization was done in [181 in case of bounded domain. 

The paper is organized as follows. In Section [5] we recall basic notions and fix notation which 
will be used throughout this paper. In Section [3] we introduce dimensionless quantities, derive ([TJ 
and (jH and set up the Cauchy problem (j4]), ([2]) in the setting of distribution theory. Section |4] 
is devoted to examination of the Cauchy problem (|4]), ([2]) where we prove a theorem on existence 
and uniqueness of its fundamental solution. In Section [5] we give a numerical example, which 
illustrates the behavior of the fundamental solution for the specific choice of the parameters a and 
T, as well as initial conditions. 



2 Mathematical preliminaries 

Let ft be an open set in M". The space of distributions will be denoted by 2?' (17) and the space 
of Schwartz tempered distributions by 5'(M"). I?'^(M) and 5:^(M) are the subspaces of 2?'(M) and 
5'(R) respectively, containing distributions supported on [0,00). Since we shall mostly work with 
functions and distributions depending on two variables, u = u(x,t), we introduce iS'(M x R^) to 
be the space of all distributions u £ 5'(R^), which vanish for t < 0. 
We shall also need the set C+ = {z e C | Re z > 0}. 

Let t € [0,a], a > 0, and y G AC([0, a]). The left Riemann-Liouville fractional derivative of 
order a £ [0, 1), oDfy, is defined as 

where F is the Euler gamma function. 

In the distributional setting, one introduces a family {fa] am. G ^'^{R) as 

f(t)^[ ^(')S' " > 0' 

1 ^rfa+Nit), a<0,a + N >0,N £N ' 

where H is the Heaviside function. Then fa* is a convolution operator acting on 2?^(R) (also, 
fa* ■ S'^{W) — > 5!^(R)). For a < it is called the operator of fractional differentiation. Moreover, 
for y £ AC{[0,a]) it coincides with the left Riemann-Liouville fractional derivative, i.e. 

oDtV ^ f-a * y- 
For y £ 5'(R) the Fourier transform is defined as 

where for G S{M) 

poo 

J^'fiO - 'fiO = / </?(a;)e-*«^ dx, f e R. 



Let y £ D'_^{M.) such that e £ S'(R), for all ^ > a > 0. Then the Laplace transform of y is 
defined by 

^y{s) ^y{s) = J^{e^^*y){v), .s = ^ + iry. 

It is well known that the function Cy is holomorphic in the half plane Re s > a (see e.g. [S] or 
pi]). In particular, for y £ L^{R) such that y{t) = 0, for t < 0, and \y{t)\ < Ae"* {a, A > 0) the 
Laplace transform is 



/•oo 

Cy{s)^ / y{t)e^'*dt, Res > 0. 
Jo 
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We shall need the following properties: y,yi,y2 € 



J" 



-y 



(0 = 



C[oDfy]{s) = 

'C[yi*y2](s) = /:yi(s) • £?;i(s), 
£(5(3) = 1. 

In order to introduce the inverse Laplace transform we recall from [21J the following: Let Y be 
a holomorphic function in the half plane Res > a such that |i^(s)| < ^■^J^J^TJiTi m,k € M. Then 
there exists a distribution y S 5!^(R) such that £y = Y, and 

y{t)^C-^Y{t) = — Y{s)e''ds, t > 0. 

J a-ioo 

The notion of fundamental solution is introduced as follows (see e.g. [21] )• Let P be a linear 
partial integro-differential operator with constant coefficients. A fundamental solution of P, de- 
noted by is a solution to the equation Pu = 6. Once the fundamental solution is determined 
one finds a solution to Pu = f a.s u = E*f. 

The Cauchy problem for the second order linear partial integro-differential operator with con- 
stant coefficients P is given by 

Pu{x,t) ^ f{x,t) (9) 
d 

u{x,0) = uo{x), —u{x,0)^ui{x), (10) 

where / is continuous for t > 0, uq G C^(M) and ui € C(]R). A classical solution u{x,t) to the 
Cauchy problem (|91IT0| is of class for t > 0, of class for t > 0, satisfies equation ^ for 
t > 0, and initial conditions (jlOp when i — !■ 0. If functions u and / are continued by zero for t < 0, 
then the following equation is satisfied in 2?'(]R.^): 

Pu = f{x,t) + uoix)S' (t) + ui{x)S{t). (11) 

The classical solutions of the Cauchy problem (|9]fT0|) are among those solutions of equation 
(fTTjl that vanish for t < 0. Therefore, the problem of finding generalized solutions (in 2?'(R^)) of 
the equation (|11|) that vanish for t < will be called generalized Cauchy problem for the operator 
P. If there is a fundamental solution E of the operator P and if / S D'{M.'^) vanishes for t < 
then there exists a unique solution to corresponding generalized Cauchy problem and is given by 

u^ E* {f{x, t) + ua{x)5'{t) ^ ui{x)5{t)). 

Again, we refer to [S], [20] or [21] for more details. 



3 Set up of the Cauchy problem 

In order to obtain ([1]), we introduce dimensionless coordinates in the system ([MZl), where ([6]) is 
modified by the fractional Zener model of the viscoelastic body: 

a{x,t)+T„fiD'^a{x,t) ^ E[e{x,t) + TefiD'^e{x,t)], xe]R,t>0, (12) 

where t„, are relaxation times satisfying > t„ > 0. The latter condition follows from the 
Second Law of Thermodynamics (see [1_ ) . 

We subject to the system (O, (O and initial conditions 

d 

u{x,{)) ^ Uf){x), —u{x,0)^vq{x), 
cr(a;,0)=0, £(a;,0)=0. 
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where x e M, uq and vq are initial displacement and velocity. Note that there are no initial stress 
and strain. We also supply boundary conditions 

lim u{x,t)—0^ lim a{x,t)=0^ t>0. 

We will need the following lemma. 

Lemma 3.1. Let y £ AC{[0,a]), a>0,0<a<l and T* > 0. Let t and y be transformed 
as {t,y) = (^,2/), then the left Riemann-Liouville fractional derivative oD"y{t), t > 0, is trans- 
formed as follows 

oDfyit) = iT*)%Dfyit). 
Proof. Since y{t) = y{t) — y{tT*), we have 



r(i - a) dt Jo - 

(TToD^yit). 



' Til^a)dtJo {t~z)-' 



□ 



Let < a < 1, > Ta > and E be constants appearing in ([T2|) . and let p be the constant in 
(O. Set T* — ^/t^ and X* — ^fr^^J^. Then X* and T* are constants (measured in meters and 

seconds respectively) which satisfy (^^^ "B = ^ ^^'^ _ 
We introduce dimensionless quantities in ([5]), ([T^ and (0 as 

_a; - t _ u _ a 

^"x^' *"t^' """x^' '^^E 

and 

Wo - - To- 

-0 = ^; -o = -^T; r = - 

and obtain system ([T][3]). Note that e is already dimensionless quantity. 

In order to simplify notation, a bar ~ will be dropped and dimensionless quantities: x, t, a, 
u, f, Mo and vq will be written as: x, t, a, u, r, uq and vq. Obviously, condition < < implies 

< T < 1. 

In our work an appropriate setting for studying system ([TJ2]) will be distributional one. In 
fact, we will look for a fundamental solution to the generalized Cauchy problem for FZWE ^ in 
S'{M. X M+). This suffices to obtain a solution to ([T][3]), since ([1]) and (|4]) are equivalent. Indeed, 
by applying the Laplace transform with respect to time variable t to the second equation in ([T]), 
one obtains 

(1 + Ts")^{x, s) = (1 + s")e(x, s). 
According to [TS], ^ i+tb" ) well-defined element in 5^(M), hence 

<^ = ^-'(^^)^tS. (13) 



1 -h TS"_ 

Inserting e from the third equation in ^ into (|13p and then a into the first equation in ([l} yields 
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Setting L{t) = C ^ ( i+ts" ) ** come to ([4]). Therefore we have proved that ([T]) and (|4]) are 
equivalent. Notice that equation ([Ti)) is of the form = 0, with 

P-S,-C-^(^):^.- (15, 



4 The existence and uniqueness of a solution to the Cauchy 
problem (jl), 

The aim of this section is to find a solution to the generalized Cauchy problem to ([4]), i.e. 

J2 / 1 _L oQ \ «2 



or equivalently 

Pu(a;, t) = uo{x)6'{t) + VQ{x)S{t), 

where P is given by (fT5|) . 

Remark 4.1. We have already explained that initial conditions are included into the generalized 
Cauchy problem (see Section (2]). For functions in S' boundary conditions @ are automatically 
fulfilled. 

We state the main theorem. 

Theorem 4.2. Let ua,vo e 5'(R). Then there exists a unique solution u £ S'{M. x M+) to U6\} 
given by 

uix, t) = S(x, t) {uQ{x)5'{t) + voix)6{t)), (17) 

where 

1 



S{x,t) = 1 + 





/l + Tq°'e^"'^ 






1 + 





(18) 



is t/ie fundamental solution of operator P, S £ 5'(IR x mi/i support in the cone \x\ < 

We will need the following lemma. 
Lemma 4.3. Let f G 5'(M). Then the equation 

v" ~ ujv = -f (19) 
has a solution v G 5'(M) for all io £ C \ (— oo, 0], which is of the form 

g-Vw|2;| 

where ^/cj is the main branch. 

Proof. We first apply Fourier transform to (|T9)) : 

(20) 
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In order to find v we need to calculate 

Let a; > 0. Let r+ = FiUFa, where Fi = {z = Re^'^ | i? > 0, < < tt} and F2 = {z = CI < 
C < R}. By the Cauchy residual theorem we have that 



/ ^txz__}_ — _ 27rzy^ Res(zj), 



where Zi are poles of the function z i-> -^2^ lying inside of F. For e C \ (—00, 0] there is only 
one such pole zi = iy/uj and Res(zi) = = '^2iJ^ • Letting R — > 00, the integral over Fi 



tends to zero and L e^^^^. — dz = L e"^ ci\ df, thus 



+ujj 27r z2 + cj 2^ 

For X < and F_ F3 U Teq-A, where F3 = {z = i?e^'^ | i? > 0, -tt < < 0} and F4 = 
= — CI ~ R < £, < R} we can apply the same arguments as for the case x > to obtain 



2yZJ ■ 

Thus for X e M we can write 



: R we can write 

Now, the inverse Fourier transform applied to (|20l) yields the result. □ 



^-^/uj\x\ 



Proof of the main theorem. Applying Laplace transform with respect to t to (|T6|) we 
obtain: 

g^uix, s) - s -j-—^u{x, s) = — ^-^—^{suo(x) + vo{x)). (21) 

Set aj(s) := \"!|T/c and/(a;,s) := (^^0(3^) +^o(a^))) s e C+. In order to applv Lemma 1431 to 

pT|) we have to show first that /(•, s) e 5'(R), which follows from assumptions that uo,vq £ S'{M.) 
and second, that uj{s) S C \ (— oo,0], for all s E C+. The latter can be verified in the following 
way. Take arbitrary s = pe"^, p>0, -§<</?< f. Then uj{s) = p2g2iy ^+^p^^^^ ^ and after a 
straightforward calculation we arrive to 



Rew(s) = ^ [(1 + p"(l + t) cos(acp) + V") cos(2(^) + p"(l - t) sin(a(^) sin(2(p)] , (22) 

Ima;(s) = ^ [(l + p"(l + r) cos(a(y5) + rp^") sm{2ip) - p"(l - r) sin(a(^) cos(2v3)] , (23) 

where A = {1 + p" cos{aip)y + p^" sm'^{aip) > 0. Suppose that uj{s) e (— cx3, 0], for some s £ C+. 
Then Ima;(s) = 0, and ^ yields 

(1 + p"(l + r) cos(a(^) + rp2") = p"(l - r) sin(a^)^^§T- (24) 
^ ' sm(2(y9) 

From (f22|) we now obtain 

Re^(s) = ^p"(l _ r) sin(a(p) sin(2^) f ^^^^ + l) < 0. 
A \sm (2(pj J 
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where we have used (|24p and the assumption uj{s) e (—00, 0], for all s e C+. 

Since - r) sin(a^) (^^^^ + l) > 0, for p > 0, -f < < f , it follows that sin(2^) < 

0, and therefore ip € (— f , — j] U [-, However, (|24|) can not be satisfied for those f. (Indeed, 
for (f G (fi~f]j cos(a(^), cos(2(/5), sin(aiy9) > 0, but sm{2ip) < and similarly for ip e [f;f)-) 
Hence, a;(s) e C \ (—00, 0], for all s G C4.. 

Now we can apply Lemma [4.31 to obtain a solution of (|2l]) : 



u{x, s) 



4s)\x\ 



uj{s) 



{suo{x) + voix)) 



2s2 



¥j: {sUo{x) + Vo{x)). 



(25) 



To obtain a solution u to (jl]) we need to calculate inverse Laplace transform of For that 

purpose set 



S{x, s) 



j{s)\x\ 



1 + 



and 



2s2 2s V 1 + 

S{x,t)^£-^ S{x,s) {t), xeM.,t>0. 
.(i+lfT) 



X eK,s e 



(26) 



S{x,t) is well-defined since |iS'(a;, s)| < ^j^^' ^ , m e N, so its inverse Laplace transform exists 
(cf. Section ID). 

Note also that S" is a fundamental solution of P. Then the solution u is given by (|17p . 
Multiform function S, given by (|26p. has branch points at s = and s = 00 and has no 
singularities. Hence, it can be evaluated by the use of the Cauchy integral formula 



S'(a;, s)e"* ds = 0, xeR,t>0, 
where F = Fi U r2 U U Fa U r4 U 70, is a contour given in Figure [1] 

tlm s 



(27) 




Figure 1: Integration contour F 
More precisely, for arbitrarily chosen i?>0,0<e<i? and a > 0, F is defined by 



Fi 


s 


= Re"^,ipo < ip <Tr; 


F2 


s 


= ge", -R < -q < -e; 


re 


s 


— ee**^, — TT < p < tt; 


r3 


s 


= ge~",e < q < R; 


F4 


s 


= Re''^, -TT < ip < -ipo; 


70 


s 


= a(l + i tan p), —ipo < p < ipo 
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where ipo = arccos(-^). Note that hmjj„j.oo i^sq = f- In the hmit when R — > oo, integral along 
contour Fi reads (a; e M, t > 0) 



i?,~i.oo 



lim / S{x, s)e'^* ds — — lini 



2 R- 



— : e V i+« = I dip. (28) 



lim 

-R— foo 



1 + ri?"e*"'^ 



-\x\Re 



Evaluating the absolute value of /p^ S{x, s)e^* ds one obtains G M, i > 0) 
Six, s)e'* ds <- lim / 

2 fl-foo 

In the limit when R oo the expression w tends to ^/r and therefore 



1 + R°'e^'^v 



1 + -R°c»°'^' 



COS 



lim 



S[x,s)e'^ ds 



< ^/f lim / e^™'''^(*-l="l^)ci(^ = 0, ifi>|a;|V^, 
2 «^o%o 



since cos(/? < for G (f '"")■ Similar argument is valid for the integral along 

In the limit when £ — > 0, the integral along is given by the formula similar to and it is 
calculated as 



1 



lim / Six, s)e ds — - lim 



1 + Te°'e^°"f -\x\ee' 

1 + £"e'"^ ^ 



i = — ZTT. 



Integrals along contours: F2, F3 and 70, in the limit when R — > 00, e 0, give 



lim / S(x, s)e ds 



lim / S{x, s)e'** ds = 



1 + Tq°'e"^-" ^-'q[t-\x\^i± 



2 7o V 1 + 9"e*"^ 



e-yO 



Fa 



1 



1 + q°'e- 



lim / 5(2;, 3)6"* ds = 27riS'(a;,t). 

Now, by the Cauchy integral formula ^7} . we obtain S as in ([T5)) . Therefore we have proved the 
theorem. □ 

As a consequence of what we proved in Sections |3] and SI we have the following corollary. 

Corollary 4.4. Let u he given by |J7p . Then 

{u,e,a){x,t)^ (^u{x,t),-^u{x,t),C-^ (l^^) ** ^"^'''^^) ^ ('^'(K >< 

is a unique solution to the system 

In a similar way one can also consider nonhomogeneous case. 
Remark 4.5. Let the first equation in ([T]) be replaced by 

d 

— cr(a;, t) Q^uix, t) + /(a;, t), 

where / G 5'(K x K+). This is a case of a rod under the influence of body forces. Then, the 
solution of the generalized Cauchy problem 



dt^ 



/ 1 + s" \ d"^ 

u{x, t) = \———\ *i -Q^uix, t) + fix, t) + uo{x)S'{t) + vo{x)5{t). 
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is 

u{x, t) = S{x, t) (/(a;, t) + ua{x)5' [t] + v^{x)5{t)) 

where S is as in (IT51) . 

Remark 4.6. Dimensionless condition < or |a;| < \f^^^i in the dimensional form, can 

be physically interpreted as the wave property. Namely, in the moment wave caused by the 
initial disturbance uq{x) = b(x)^ a; € R, has reached the point which is at distance |x| from the 

origin, where the initial disturbance was applied. Thus, constant cz — \J~E\j^ ^^"^ interpreted 

as the wave speed in the fractional viscoelastic media of Zener type. 



5 A numerical example 

Let uq — b and f o = in the solution (jl7p to the Cauchy problem of wave equation for fractional 
Zener type viscoelastic media. Then solution reads 

d 

u{x, t) = -qiS{x, t), a; e M, <> 0, (29) 

i.e. the fundamental solution represents the solution itself. Figure [5] presents the plots of u{x,t), 
X e (0, 3), t e {0.5, 1, 1.5}, given by for the set of parameters: a = 0.23 and t = 0.004. 

u(x,t) 




0.5 1 1.5 2 2.5 3 



Figure 2: Solution u{x, t), x e{0,3),t£ {0.5, 1, 1.5} 



In order to see the effect of a change in the order of a fractional derivative, we present following 
figures. Figure |3] presents plots oi u{x,t), x e (0,3), t e {0.5,1,1.5}, given by (|29|) . for different 
values of parameter a. Plot of u, denoted by the dashed line corresponds to a = 0.25, dot-dashed 
line is used for a — 0.5, while the solid line denotes the plot for a = 0.75. Parameter r = 0.004 is 
the same in all three figures. 



U(x,t) ri(x,t) u(x,t) 




Figure 3: Solution u{x, t), x £ (0, 6), t £ {0.5, 1, 1.5}, a £ {0.25, 0.5, 0.75} 

From figure [3] one can see that for each value of the order of fractional derivative a, as time 
increases, the maximum value of u decreases, which is the consequence of the dissipative model of 
a media. 
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Figure m presents plots of u{x,t), x e (0,3), a € {0.25,0.5,0.75}, given by (HH), for different 
time moments, i.e. for t G {0.5, 1, 1.5}. Parameter r and line styles are the same as in figure|31 



u(x,0.5) u(x,l) u{x,1.5) 




Figure 4: Solution t), x e (0, 5), t e {0.5, 1, 1.5}, a e {0.25, 0.5, 0.75} 

From figure |4] one can see that at the same time instance, as the value of a increases, the 
maximum value of u decreases, while the width of the maximum increases. This is the consequence 
of the fact that the constitutive equation of the fractional type in ([T]) describes media that tends 
to be elastic as a tends to zero. Therefore, the fundamental solution ((29)) should tend to the Dirac 
distribution as a — 0. Also, media tends to be viscoelastic as a in ([T|) tends to 1. Therefore, the 
fundamental solution ([29]) should tend to the solution of wave equation for viscoelastic media as 
a — > 1. 

Figure [S] presents plots of u for both different time moments and values of a. Parameter r and 
line styles are the same as in figure [3l 

u(x,t) 

1.2[ /'. t = 0.5 

1 

0.8 
0.6 

. 4 
. 2 



Figure 5: Solution u{x, t), x G (0, 5), t G {0.5, 1, 1.5}, a G {0.25, 0.5, 0.75} 
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